function [out, zlb_f, z_f_T, Q_out, G_out] = zlb(SET, in, e_in, Qt_in, Gt_in, FG)
% ZLB as piecewise-linear approximation
%
% jonescallum@gmail.com

if nargin<4
    Qt_in = repmat(SET.mats.Q,[1 1 SET.maxchk]) ;
    Gt_in = repmat(SET.mats.G,[1 1 SET.maxchk]) ;
end

if nargin<6
    z_f_i = 100 ;
    z_f_f = 0 ;
end

if nargin==6
    z_f_i = FG.z_f_i ; 
    z_f_f = FG.z_f_f ; 
end

%%

i_f     = SET.pos_of_i ;
zlb_val = SET.zlb_val ;
maxchk  = SET.maxchk ;

Qt = Qt_in;%zeros(n_,n_,1);
Gt = Gt_in;%zeros(n_,l_,1);

%% Step 1(a)

ychk = Qt(:,:,1) * in + Gt(:,:,1) * e_in ;

for t_ = 2:maxchk
    ychk(:,t_) = Qt(:,:,t_) * ychk(:,t_-1) ;
end

%% Step 1(b)

chg_f = ychk(i_f,:)<zlb_val-1e-10 ;
chk_f = sum(chg_f) ;
chk   = chk_f ;

if ~chk
    switch SET.stoch
        case 1 ; out = ychk(:,1) ;
        case 0 ; out = ychk(:,1:maxchk) ;
    end
    zlb_f = 0 ; 
    z_f_T = 0 ; 
    Q_out = Qt(:,:,1) ;
    G_out = Gt(:,:,1) ;    
    return
end

%% Step 1(c)

max_t = z_f_f ;

zlb_f = chk_f>0 ; 

while chk_f

    for t_ = 1:length(chg_f)
        if chg_f(t_)==1
            
            z_f_i = min(z_f_i,t_) ;
            z_f_f = max(z_f_f,t_) ;
            max_t = max(max_t, z_f_f) ;
            
            SET.tv_str_chg.Qf = Qt_in(:,:,z_f_f+1) ;
            
            mats = rfmats(SET, z_f_i, z_f_f) ;
            
            Qt_ = mats.Qhat ;
            Gt_ = mats.Ghat ;
            
            break
        end
    end
    
    ychk(:,1) = Qt_(:,:,1) * in + Gt_(:,:,1) * e_in ;
    
    for tt = 2:z_f_f
        ychk(:,tt) = Qt_(:,:,tt) * ychk(:,tt-1) ;
    end
        
    for tt = (z_f_f+1):maxchk
        ychk(:,tt) = Qt(:,:,tt) * ychk(:,tt-1) ;
    end
    
chg_f = ychk(i_f,:)<zlb_val-1e-10 ;    

chk_f = sum(chg_f)>0 ;

end

if zlb_f ; z_f_T = z_f_f - z_f_i + 1 ;
else       z_f_T = 0 ;
end

%% Outputs

switch SET.stoch
    case 1 ; out = ychk(:,1) ; 
    case 0 ; out = ychk(:,1:maxchk) ;     
end

Q_out = Qt_(:,:,1) ;
G_out = Gt_(:,:,1) ;